# GPS L1Cp code construction
#
# Copyright 2018 Peter Monta

import numpy as np
from sympy.ntheory import legendre_symbol

chip_rate = 1023000
code_length = 10230

l1cp_params = {
   1: (5111,412),    2: (5109,161),     3: (5108,1),      4: (5106,303),
   5: (5103,207),    6: (5101,4971),    7: (5100,4496),   8: (5098,5),
   9: (5095,4557),  10: (5094,485),    11: (5093,253),   12: (5091,4676),
  13: (5090,1),     14: (5081,66),     15: (5080,4485),  16: (5069,282),
  17: (5068,193),   18: (5054,5211),   19: (5044,729),   20: (5027,4848),
  21: (5026,982),   22: (5014,5955),   23: (5004,9805),  24: (4980,670),
  25: (4915,464),   26: (4909,29),     27: (4893,429),   28: (4885,394),
  29: (4832,616),   30: (4824,9457),   31: (4591,4429),  32: (3706,4771),
  33: (5092,365),   34: (4986,9705),   35: (4965,9489),  36: (4920,4193),
  37: (4917,9947),  38: (4858,824),    39: (4847,864),   40: (4790,347),
  41: (4770,677),   42: (4318,6544),   43: (4126,6312),  44: (3961,9804),
  45: (3790,278),   46: (4911,9461),   47: (4881,444),   48: (4827,4839),
  49: (4795,4144),  50: (4789,9875),   51: (4725,197),   52: (4675,1156),
  53: (4539,4674),  54: (4535,10035),  55: (4458,4504),  56: (4197,5),
  57: (4096,9937),  58: (3484,430),    59: (3481,5),     60: (3393,355),
  61: (3175,909),   62: (2360,1622),   63: (1852,6284),
  64: (5065,9429),  65: (5063,77),     66: (5055,932),   67: (5012,5973),
  68: (4981,377),   69: (4952,10000),  70: (4934,951),   71: (4932,6212),
  72: (4786,686),   73: (4762,9352),   74: (4640,5999),  75: (4601,9912),
  76: (4563,9620),  77: (4388,635),    78: (3820,4951),  79: (3687,5453),
  80: (5052,4658),  81: (5051,4800),   82: (5047,59),    83: (5039,318),
  84: (5015,571),   85: (5005,565),    86: (4984,9947),  87: (4975,4654),
  88: (4974,148),   89: (4972,3929),   90: (4962,293),   91: (4913,178),
  92: (4907,10142), 93: (4903,9683),   94: (4833,137),   95: (4778,565),
  96: (4721,35),    97: (4661,5949),   98: (4660,2),     99: (4655,5982),
 100: (4623,825),  101: (4590,9614),  102: (4548,9790), 103: (4461,5613),
 104: (4442,764),  105: (4347,660),   106: (4259,4870), 107: (4256,4950),
 108: (4166,4881), 109: (4155,1151),  110: (4109,9977), 111: (4100,5122),
 112: (4023,10074),113: (3998,4832),  114: (3979,77),   115: (3903,4698),
 116: (3568,1002), 117: (5088,5549),  118: (5050,9606), 119: (5020,9228),
 120: (4990,604),  121: (4982,4678),  122: (4966,4854), 123: (4949,4122),
 124: (4947,9471), 125: (4937,5026),  126: (4935,272),  127: (4906,1027),
 128: (4901,317),  129: (4872,691),   130: (4865,509),  131: (4863,9708),
 132: (4818,5033), 133: (4785,9938),  134: (4781,4314), 135: (4776,10140),
 136: (4775,4790), 137: (4754,9823),  138: (4696,6093), 139: (4690,469),
 140: (4658,1215), 141: (4607,799),   142: (4599,756),  143: (4596,9994),
 144: (4530,4843), 145: (4524,5271),  146: (4451,9661), 147: (4441,6255),
 148: (4396,5203), 149: (4340,203),   150: (4335,10070),151: (4296,30),
 152: (4267,103),  153: (4168,5692),  154: (4149,32),   155: (4097,9826),
 156: (4061,76),   157: (3989,59),    158: (3966,6831), 159: (3789,958),
 160: (3775,1471), 161: (3622,10070), 162: (3523,553),  163: (3515,5487),
 164: (3492,55),   165: (3345,208),   166: (3235,645),  167: (3169,5268),
 168: (3157,1873), 169: (3082,427),   170: (3072,367),  171: (3032,1404),
 172: (3030,5652), 173: (4582,5),     174: (4595,368),  175: (4068,451),
 176: (4871,9595), 177: (4514,1030),  178: (4439,1324), 179: (4122,692),
 180: (4948,9819), 181: (4774,4520),  182: (3923,9911), 183: (3411,278),
 184: (4745,642),  185: (4195,6330),  186: (4897,5508), 187: (3047,1872),
 188: (4185,5445), 189: (4354,10131), 190: (5077,422),  191: (4042,4918),
 192: (2111,787),  193: (4311,9864),  194: (5024,9753), 195: (4352,9859),
 196: (4678,328),  197: (5034,1),     198: (5085,4733), 199: (3646,164),
 200: (4868,135),  201: (3668,174),   202: (4211,132),  203: (2883,538),
 204: (2850,176),  205: (2815,198),   206: (2542,595),  207: (2492,574),
 208: (2376,321),  209: (2036,596),   210: (1920,491),
}

N = 10223
L = np.array([legendre_symbol(i,N) for i in range(N)])
L[L==-1] = 0
L[0] = 0

def l1cp(prn):
  w,p = l1cp_params[prn]
  W = np.array([L[k]^L[(k+w)%N] for k in range(N)])
  expansion = np.array([0,1,1,0,1,0,0])
  c = np.concatenate((W[0:p-1],expansion,W[p-1:N]))
  return c

codes = {}

def l1cp_code(prn):
  if prn not in codes:
    codes[prn] = l1cp(prn)
  return codes[prn]

def code(prn,chips,frac,incr,n):
  c = l1cp_code(prn)
  idx = (chips%code_length) + frac + incr*np.arange(n)
  idx = np.floor(idx).astype('int')
  idx = np.mod(idx,code_length)
  x = c[idx]
  return 1.0 - 2.0*x

l1cp_secondary_params = {
   1: (0o5111,0o3266),   2: (0o5421,0o2040),   3: (0o5501,0o1527),   4: (0o5403,0o3307),
   5: (0o6417,0o3756),   6: (0o6141,0o3026),   7: (0o6351,0o0562),   8: (0o6501,0o0420),
   9: (0o6205,0o3415),  10: (0o6235,0o0337),  11: (0o7751,0o0265),  12: (0o6623,0o1230),
  13: (0o6733,0o2204),  14: (0o7627,0o1440),  15: (0o5667,0o2412),  16: (0o5051,0o3516),
  17: (0o7665,0o2761),  18: (0o6325,0o3750),  19: (0o4365,0o2701),  20: (0o4745,0o1206),
  21: (0o7633,0o1544),  22: (0o6747,0o1774),  23: (0o4475,0o0546),  24: (0o4225,0o2213),
  25: (0o7063,0o3707),  26: (0o4423,0o2051),  27: (0o6651,0o3650),  28: (0o4161,0o1777),
  29: (0o7237,0o3203),  30: (0o4473,0o1762),  31: (0o5477,0o2100),  32: (0o6163,0o0571),
  33: (0o7223,0o3710),  34: (0o6323,0o3535),  35: (0o7125,0o3110),  36: (0o7035,0o1426),
  37: (0o4341,0o0255),  38: (0o4353,0o0321),  39: (0o4107,0o3124),  40: (0o5735,0o0572),
  41: (0o6741,0o1736),  42: (0o7071,0o3306),  43: (0o4563,0o1307),  44: (0o5755,0o3763),
  45: (0o6127,0o1604),  46: (0o4671,0o1021),  47: (0o4511,0o2624),  48: (0o4533,0o0406),
  49: (0o5357,0o0114),  50: (0o5607,0o0077),  51: (0o6673,0o3477),  52: (0o6153,0o1000),
  53: (0o7565,0o3460),  54: (0o7107,0o2607),  55: (0o6211,0o2057),  56: (0o4321,0o3467),
  57: (0o7201,0o0706),  58: (0o4451,0o2032),  59: (0o5411,0o1464),  60: (0o5141,0o0520),
  61: (0o7041,0o1766),  62: (0o6637,0o3270),  63: (0o4577,0o0341),
  64: (0o5111,0o1740,0o3035),  65: (0o5111,0o3664,0o1557),  66: (0o5111,0o1427,0o0237),  67: (0o5111,0o2627,0o2527),
  68: (0o5111,0o0701,0o3307),  69: (0o5111,0o3460,0o1402),  70: (0o5111,0o1373,0o1225),  71: (0o5111,0o2540,0o0607),
  72: (0o5111,0o2004,0o0351),  73: (0o5111,0o2274,0o3724),  74: (0o5111,0o1340,0o1675),  75: (0o5111,0o0602,0o2625),
  76: (0o5111,0o2502,0o1030),  77: (0o5111,0o0327,0o1443),  78: (0o5111,0o2600,0o3277),  79: (0o5111,0o0464,0o1132),
  80: (0o5111,0o3674,0o0572),  81: (0o5111,0o3040,0o1241),  82: (0o5111,0o1153,0o0535),  83: (0o5111,0o0747,0o1366),
  84: (0o5111,0o1770,0o0041),  85: (0o5111,0o3772,0o0561),  86: (0o5111,0o1731,0o0122),  87: (0o5111,0o1672,0o1205),
  88: (0o5111,0o1333,0o3753),  89: (0o5111,0o2705,0o2543),  90: (0o5111,0o2713,0o3031),  91: (0o5111,0o3562,0o2260),
  92: (0o5111,0o3245,0o3773),  93: (0o5111,0o3770,0o3156),  94: (0o5111,0o3202,0o2215),  95: (0o5111,0o3521,0o0146),
  96: (0o5111,0o3250,0o2413),  97: (0o5111,0o2117,0o2564),  98: (0o5111,0o0530,0o3310),  99: (0o5111,0o3021,0o2267),
 100: (0o5421,0o2511,0o3120), 101: (0o5421,0o1562,0o0064), 102: (0o5421,0o1067,0o1042), 103: (0o5421,0o0424,0o0476),
 104: (0o5421,0o3402,0o1020), 105: (0o5421,0o1326,0o0431), 106: (0o5421,0o2142,0o0216), 107: (0o5421,0o0733,0o2736),
 108: (0o5421,0o0504,0o2527), 109: (0o5421,0o1611,0o2431), 110: (0o5421,0o2724,0o1013), 111: (0o5421,0o0753,0o0524),
 112: (0o5421,0o3724,0o0726), 113: (0o5421,0o2652,0o1042), 114: (0o5421,0o1743,0o3362), 115: (0o5421,0o0013,0o1364),
 116: (0o5421,0o3464,0o3354), 117: (0o5421,0o2300,0o0623), 118: (0o5421,0o1334,0o0145), 119: (0o5421,0o2175,0o0214),
 120: (0o5421,0o2564,0o0223), 121: (0o5421,0o3075,0o0151), 122: (0o5421,0o3455,0o2405), 123: (0o5421,0o3627,0o2522),
 124: (0o5421,0o0617,0o3235), 125: (0o5421,0o1324,0o0452), 126: (0o5421,0o3506,0o2617), 127: (0o5421,0o2231,0o1300),
 128: (0o5421,0o1110,0o1430), 129: (0o5421,0o1271,0o0773), 130: (0o5421,0o3740,0o0772), 131: (0o5421,0o3652,0o3561),
 132: (0o5421,0o1644,0o0607), 133: (0o5421,0o3635,0o0420), 134: (0o5421,0o3436,0o0527), 135: (0o5421,0o3076,0o3770),
 136: (0o5421,0o0434,0o2536), 137: (0o5421,0o3340,0o2233), 138: (0o5421,0o0054,0o3366), 139: (0o5403,0o2446,0o3766),
 140: (0o5403,0o0025,0o3554), 141: (0o5403,0o0150,0o2060), 142: (0o5403,0o2746,0o2070), 143: (0o5403,0o2723,0o0713),
 144: (0o5403,0o2601,0o3366), 145: (0o5403,0o3440,0o3247), 146: (0o5403,0o1312,0o2776), 147: (0o5403,0o0544,0o1244),
 148: (0o5403,0o2062,0o2102), 149: (0o5403,0o0176,0o1712), 150: (0o5403,0o3616,0o1245), 151: (0o5403,0o1740,0o3344),
 152: (0o5403,0o3777,0o1277), 153: (0o5403,0o0432,0o0165), 154: (0o5403,0o2466,0o2131), 155: (0o5403,0o1667,0o3623),
 156: (0o5403,0o3601,0o0141), 157: (0o5403,0o2706,0o0421), 158: (0o5403,0o2022,0o3032), 159: (0o5403,0o1363,0o2065),
 160: (0o5403,0o2331,0o3024), 161: (0o5403,0o3556,0o2663), 162: (0o5403,0o2205,0o2274), 163: (0o5403,0o3734,0o2114),
 164: (0o5403,0o2115,0o1664), 165: (0o5403,0o0010,0o0413), 166: (0o5403,0o2140,0o1512), 167: (0o5403,0o3136,0o0135),
 168: (0o5403,0o0272,0o2737), 169: (0o5403,0o3264,0o1015), 170: (0o5403,0o2017,0o1075), 171: (0o5403,0o2505,0o1255),
 172: (0o5403,0o3532,0o3473), 173: (0o5403,0o0647,0o2716), 174: (0o5403,0o1542,0o0101), 175: (0o5403,0o2154,0o1105),
 176: (0o5403,0o3734,0o1407), 177: (0o5403,0o2621,0o3407), 178: (0o5403,0o2711,0o1046), 179: (0o5403,0o0217,0o3237),
 180: (0o5403,0o3503,0o0154), 181: (0o5403,0o3457,0o3010), 182: (0o5403,0o3750,0o2245), 183: (0o5403,0o2525,0o2051),
 184: (0o5403,0o0113,0o2144), 185: (0o5403,0o0265,0o1743), 186: (0o5403,0o1711,0o2511), 187: (0o5403,0o0552,0o3410),
 188: (0o5403,0o0675,0o1414), 189: (0o5403,0o1706,0o1275), 190: (0o5403,0o3513,0o2257), 191: (0o5403,0o1135,0o2331),
 192: (0o5403,0o0566,0o0276), 193: (0o5403,0o0500,0o3261), 194: (0o5403,0o0254,0o1760), 195: (0o5403,0o3445,0o0430),
 196: (0o5403,0o2542,0o3477), 197: (0o5403,0o1257,0o1676), 198: (0o6501,0o0211,0o1636), 199: (0o6501,0o0534,0o2411),
 200: (0o6501,0o1420,0o1473), 201: (0o6501,0o3401,0o2266), 202: (0o6501,0o0714,0o2104), 203: (0o6501,0o0613,0o2070),
 204: (0o6501,0o2475,0o1766), 205: (0o6501,0o2572,0o0711), 206: (0o6501,0o3265,0o2533), 207: (0o6501,0o1250,0o0353),
 208: (0o6501,0o1711,0o1744), 209: (0o6501,0o2704,0o0053), 210: (0o6501,0o0135,0o2222),
}

sec_code_length = 1800

def int2list(x,n):
  y = []
  for i in range(n):
    y.append((x>>i)&1)
  return y

def xorprod(a,b):
  t = 0
  for x,y in zip(a,b):
    t = t ^ (x*y)
  return t

def s_shift(x,p):
  return [xorprod(x,p)] + x[0:-1]

def sec_l1cp(prn):
  p,init = l1cp_secondary_params[prn]
  p = int2list(p//2,11)
  x = int2list(init,11)
  c = np.zeros(sec_code_length)
  for i in range(sec_code_length):
    c[i] = x[10]
    x = s_shift(x,p)
  return c

def sec_l1cp_extended(prn):
  p1,init1,init2 =  l1cp_secondary_params[prn]
  p2 = 0o5001
  p1 = int2list(p1//2,11)
  x1 = int2list(init1,11)
  p2 = int2list(p2//2,11)
  x2 = int2list(init2,11)
  c = np.zeros(sec_code_length)
  for i in range(sec_code_length):
    c[i] = x1[10] ^ x2[10]
    x1 = s_shift(x1,p1)
    x2 = s_shift(x2,p2)
  return c

secondary_codes = {}

def secondary_code(prn):
  if prn not in secondary_codes:
    if prn<64:
      secondary_codes[prn] = sec_l1cp(prn)
    else:
      secondary_codes[prn] = sec_l1cp_extended(prn)
  return secondary_codes[prn]

boc11 = np.array([1.0,-1.0])
tmboc_pattern = np.array([1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0])

try:
  from numba import jit
except:
  def jit(**kwargs):
    return lambda x: x

@jit(nopython=True)
def correlate(x,prn,chips,frac,incr,c,boc11):
  n = len(x)
  p = 0.0j
  cp = (chips+frac)%code_length
  bp = (2*(chips+frac))%2
  bp6 = (12*(chips+frac))%2
  u = int(cp%33)
  for i in range(n):
    if tmboc_pattern[u]:
      boc = boc11[int(bp6)]
    else:
      boc = boc11[int(bp)]
    p += x[i]*(1.0-2.0*c[int(cp)])*boc
    cp = (cp+incr)%code_length
    bp = (bp+2*incr)%2
    bp6 = (bp6+12*incr)%2
    u = int(cp%33)
  return p

# test

def chips2octal(x):
  s = ''
  for i in range(len(x)//3):
    d = 4*x[3*i] + 2*x[3*i+1] + x[3*i+2]
    s = s + '%o'%int(d)
  return s

if __name__=='__main__':
  for prn in range(1,211):
    c = l1cp_code(prn)
    s1 = chips2octal(c[0:24])
    s2 = chips2octal(c[-24:])
    print("%d %s %s"%(prn,s1,s2))
  print("secondary:")
  for prn in range(1,211):
    c = secondary_code(prn)
    print('%d %s'%(prn,chips2octal(np.concatenate((np.array([0]),c[-11:])))))
